/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include "qfunc.h"
#include "prating.c"

double
qfunc (const double*                    alpha,
       unsigned int                     n_workers,
       unsigned int                     n_labels,
       const double*                    logbeta,
       unsigned int                     n_images,
       const double*                    gamma,
       const double*                    pz,
       const double*                    priorz,
       const Rating*                    ratings,
       unsigned int                     n_ratings,
       const ParameterizedDistribution* prioralpha,
       const Distribution*              priorlogbeta,
       const Distribution*              priorgamma,
       int                              symmetric)
{
  double q = 0;

  /* TODO: This first loop computes a constant as far as the minimization
   * routine is concerned.  Maybe it can be eliminated?  (Maybe why not,
   * the termination condition is a function of the relative change in q).
   */

  for (unsigned int i = 0; i < n_images; ++i)
    {
      for (unsigned int j = 0; j < n_labels; ++j)
        {
          q += pz[j + i * n_labels] * log (priorz[j]);
        }
    }

  for (unsigned int i = 0; i < n_ratings; ++i)
    {
      for (unsigned int j = 0; j < n_labels; ++j)
        {
#define alphafunc(_k)                                           \
  ({                                                            \
    (_k == j) ? 0.0 :                                           \
    (symmetric && _k < j)                                       \
        ? alpha[ratings[i].worker * n_labels * n_labels +       \
                _k * n_labels +                                 \
                j]                                              \
        : alpha[ratings[i].worker * n_labels * n_labels +       \
                j * n_labels +                                  \
                _k];                                            \
  })

          q += pz[j + ratings[i].image * n_labels] * 
               logprating (ratings[i].rating,
                           alphafunc,
                           logbeta[ratings[i].image],
                           n_labels);
        }
#undef alphafunc
    }

  for (unsigned int i = 0; i < n_workers; ++i)
    {
      for (unsigned int j = 0; j < n_labels; ++j)
        {
          for (unsigned int k = 0; k < n_labels; ++k)
            {
              if (k == j)
                {
                  q += parameterized_distribution_logpdf (prioralpha, 1.0, 0.0);
                }
              else if (symmetric && k < j)
                {
                  q += parameterized_distribution_logpdf 
                         (prioralpha,
                          gamma[k * n_labels + j],
                          alpha[i * n_labels * n_labels +
                                k * n_labels +
                                j]);
                }
              else
                {
                  q += parameterized_distribution_logpdf 
                         (prioralpha,
                          gamma[j * n_labels + k],
                          alpha[i * n_labels * n_labels +
                                j * n_labels +
                                k]);
                }
            }
        }
    }

  for (unsigned int i = 0; i < n_images; ++i)
    {
      q += distribution_logpdf (priorlogbeta, logbeta[i]);
    }

  for (unsigned int j = 0; j < n_labels; ++j)
    {
      for (unsigned int k = 0; k < n_labels; ++k)
        {
          if (k == j)
            {
              q += distribution_logpdf (priorgamma, 1.0);
            }
          else if (symmetric && k < j)
            {
              q += distribution_logpdf (priorgamma, gamma[k * n_labels + j]);
            }
          else
            {
              q += distribution_logpdf (priorgamma, gamma[j * n_labels + k]);
            }
        }
    }

  return q;
}

void
dqfunc (double*                         dalpha,
        double*                         dlogbeta,
        double*                         dgamma,
        const double*                   alpha,
        unsigned int                    n_workers,
        unsigned int                    n_labels,
        const double*                   logbeta,
        unsigned int                    n_images,
        const double*                   gamma,
        const double*                   pz,
        const Rating*                   ratings,
        unsigned int                    n_ratings,
        const ParameterizedDistribution* prioralpha,
        const Distribution*             priorlogbeta,
        const Distribution*             priorgamma,
        int                             symmetric)
{
  memset (dalpha, 0, n_workers * n_labels * n_labels * sizeof (dalpha[0]));
  memset (dlogbeta, 0, n_images * sizeof (dlogbeta[0]));
  memset (dgamma, 0, n_labels * n_labels * sizeof (dgamma[0]));

  for (unsigned int i = 0; i < n_ratings; ++i)
    {
      for (unsigned int j = 0; j < n_labels; ++j)
        {
          double tmpdalpha[n_labels];

#define alphafunc(_k)                                           \
  ({                                                            \
    (_k == j) ? 0.0 :                                           \
    (symmetric && _k < j)                                       \
        ? alpha[ratings[i].worker * n_labels * n_labels +       \
                _k * n_labels +                                 \
                j]                                              \
        : alpha[ratings[i].worker * n_labels * n_labels +       \
                j * n_labels +                                  \
                _k];                                            \
  })

          dlogbeta[ratings[i].image] += 
            pz[j + ratings[i].image * n_labels] * 
            dlogpratingdlogbeta (ratings[i].rating,
                                 alphafunc,
                                 logbeta[ratings[i].image],
                                 n_labels);

          dlogpratingdalpha (tmpdalpha,
                             ratings[i].rating,
                             alphafunc,
                             logbeta[ratings[i].image],
                             n_labels);
#undef alphafunc

          for (unsigned int k = 0; k < n_labels; ++k)
            {
              if (k == j)
                {
                  /* do nothing ... */
                }
              else if (symmetric && k < j) 
                {
                  dalpha[ratings[i].worker * n_labels * n_labels +
                         k * n_labels +
                         j] += 
                    pz[j + ratings[i].image * n_labels] * tmpdalpha[k];
                }
              else
                {
                  dalpha[ratings[i].worker * n_labels * n_labels +
                         j * n_labels +
                         k] += 
                    pz[j + ratings[i].image * n_labels] * tmpdalpha[k];
                }
            }
        }
    }

  for (unsigned int i = 0; i < n_workers; ++i)
    {
      for (unsigned int j = 0; j < n_labels; ++j)
        {
          for (unsigned int k = 0; k < n_labels; ++k)
            {
              if (k == j)
                {
                  /* do nothing ... */
                }
              else if (symmetric && k < j) 
                {
                  dalpha[i * n_labels * n_labels +
                         k * n_labels +
                         j] += 
                    parameterized_distribution_d1logpdf
                      (prioralpha,
                       gamma[k * n_labels + j],
                       alpha[i * n_labels * n_labels +
                             k * n_labels +
                             j]);

                  dgamma[k * n_labels + j] += 
                    parameterized_distribution_d0logpdf
                      (prioralpha,
                       gamma[k * n_labels + j],
                       alpha[i * n_labels * n_labels +
                             k * n_labels +
                             j]);
                }
              else
                {
                  dalpha[i * n_labels * n_labels +
                         j * n_labels +
                         k] += 
                    parameterized_distribution_d1logpdf
                      (prioralpha,
                       gamma[j * n_labels + k],
                       alpha[i * n_labels * n_labels +
                             j * n_labels +
                             k]);

                  dgamma[j * n_labels + k] +=
                    parameterized_distribution_d0logpdf
                      (prioralpha,
                       gamma[j * n_labels + k],
                       alpha[i * n_labels * n_labels +
                             j * n_labels +
                             k]);
                }
            }
        }
    }

  for (unsigned int i = 0; i < n_images; ++i)
    {
      dlogbeta[i] += distribution_dlogpdf (priorlogbeta, logbeta[i]);
    }

  for (unsigned int j = 0; j < n_labels; ++j)
    {
      for (unsigned int k = 0; k < n_labels; ++k)
        {
          if (k == j)
            {
              /* do nothing ... */
            }
          else if (symmetric && k < j)
            {
              dgamma[k * n_labels + j] +=
                distribution_dlogpdf (priorgamma, gamma[k * n_labels + j]);
            }
          else
            {
              dgamma[j * n_labels + k] +=
                distribution_dlogpdf (priorgamma, gamma[j * n_labels + k]);
            }
        }
    }
}
